{
 "metadata": {
  "name": "",
  "signature": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os.path as osp\n",
      "\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "import scipy.linalg as la\n",
      "\n",
      "import logging\n",
      "logging.getLogger().setLevel(logging.DEBUG)\n",
      "\n",
      "import openmodes\n",
      "import openmodes.basis\n",
      "from openmodes.sources import PlaneWaveSource\n",
      "from openmodes.integration import DunavantRule\n",
      "from openmodes.parts import SinglePart\n",
      "from openmodes.mesh import TriangularSurfaceMesh\n",
      "import openmodes.gmsh"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from openmodes.ipython import init_3d\n",
      "init_3d()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Load the geometry\n",
      "\n",
      "name = 'sphere'\n",
      "parameters = {'radius': 5e-3, 'mesh_tol': 2e-3}\n",
      "meshed_name = openmodes.gmsh.mesh_geometry(osp.join(openmodes.geometry_dir, name+'.geo'), parameters['mesh_tol'], parameters=parameters)\n",
      "\n",
      "raw_mesh = openmodes.gmsh.read_mesh(meshed_name)\n",
      "\n",
      "meshes = tuple(TriangularSurfaceMesh(sub_mesh) for sub_mesh in raw_mesh)\n",
      "\n",
      "parent_part = SinglePart(meshes[0])\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create the basis functions\n",
      "\n",
      "basis_class = openmodes.basis.DivRwgBasis\n",
      "basis = basis_class(parent_part.mesh)\n",
      "\n",
      "#container = openmodes.basis.BasisContainer(basis_class)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "integration_rule = DunavantRule(12)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Set the operating frequency\n",
      "\n",
      "s= 2j*np.pi*1e9"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the source vector\n",
      "\n",
      "e_inc = np.array([1, 0, 0], dtype=np.complex128)\n",
      "k_hat = np.array([0, 0, 1], dtype=np.complex128)\n",
      "pw = PlaneWaveSource(e_inc, k_hat)\n",
      "\n",
      "field = lambda r: pw.magnetic_field(s, r)\n",
      "\n",
      "V = basis.weight_function(field, integration_rule, parent_part.nodes, n_cross=False)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import openmodes.core\n",
      "\n",
      "from openmodes.core import z_mfie_faces_self\n",
      "\n",
      "\n",
      "def impedance_rwg_mfie_free_space(s, integration_rule, basis_o, nodes_o,\n",
      "                                  basis_s, nodes_s, normals, self_impedance,\n",
      "                                  tangential_form):\n",
      "    \"\"\"MFIE derived Impedance matrix for RWG or loop-star basis functions\"\"\"\n",
      "\n",
      "    transform_o, _ = basis_o.transformation_matrices\n",
      "    num_faces_o = len(basis_o.mesh.polygons)\n",
      "\n",
      "    if self_impedance:\n",
      "        # calculate self impedance\n",
      "\n",
      "        num_faces_s = num_faces_o\n",
      "        Z_faces = z_mfie_faces_self(nodes_o, basis_o.mesh.polygons,\n",
      "                                    basis_o.mesh.polygon_areas, s,\n",
      "                                    integration_rule.xi_eta,\n",
      "                                    integration_rule.weights, normals,\n",
      "                                    tangential_form)\n",
      "\n",
      "        transform_s = transform_o\n",
      "\n",
      "    else:\n",
      "        # calculate mutual impedance\n",
      "        raise NotImplementedError\n",
      "\n",
      "        num_faces_s = len(basis_s.mesh.polygons)\n",
      "\n",
      "        transform_Z_s, _ = basis_s.transformation_matrices\n",
      "\n",
      "    if np.any(np.isnan(Z_faces)):\n",
      "        raise ValueError(\"NaN returned in impedance matrix\")\n",
      "\n",
      "    Z = transform_o.dot(transform_s.dot(Z_faces.reshape(num_faces_o*3,\n",
      "                                                        num_faces_s*3,\n",
      "                                                        order='C').T).T)\n",
      "    return Z\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Fill and solve the impedance matrix\n",
      "\n",
      "nodes = parent_part.nodes\n",
      "\n",
      "Z = impedance_rwg_mfie_free_space(s, integration_rule, basis, nodes,\n",
      "                                  basis, nodes, basis.mesh.surface_normals, self_impedance=True, tangential_form=True)\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "I = la.solve(Z, V)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "compress_scalars = None\n",
      "compress_separately = False\n",
      "\n",
      "centre, current, charge = basis.interpolate_function(I*100, return_scalar=True, nodes=parent_part.nodes)\n",
      "\n",
      "#if compress_scalars:\n",
      "#    charge = compress(charge, compress_scalars)\n",
      "\n",
      "\n",
      "\n",
      "openmodes.ipython.plot_3d([parent_part], [charge], [current], [centre])\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(10, 5))\n",
      "plt.plot(abs(I), '+')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": ""
    }
   ],
   "metadata": {}
  }
 ]
}